An Old Disk That Can Still Form a Planetary System 
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From the masses of planets orbiting our Sun, and relative elemental abundances, it is esti- 
mated that at birth our Solar System required a minimum disk mass of ~0.01 M within 
~100 AU of the star 1 " 4 . The main constituent, gaseous molecular hydrogen, does not emit 
from the disk mass reservoir^, so the most common measure of the disk mass is dust ther- 
mal emission and lines of gaseous carbon monoxide 6 . Carbon monoxide emission generally 
probes the disk surface, while the conversion from dust emission to gas mass requires knowl- 
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edge of the grain properties and gas-to-dust mass ratio, which likely differ from their inter- 
stellar values 7 8 . Thus, mass estimates vary by orders of magnitude, as exemplified by the rel- 
atively old (3-10 Myr) star TW Hya^ with estimates ranging from 0.0005 to 0.06 M ™31 
Here we report the detection the fundamental rotational transition of hydrogen deuteride, 
HD, toward TW Hya. HD is a good tracer of disk gas because it follows the distribution of 
molecular hydrogen and its emission is sensitive to the total mass. The HD detection, com- 
bined with existing observations and detailed models, implies a disk mass >0.05 M , enough 
to form a planetary system like our own. 

Commonly used tracers of protoplanetary disk masses are thermal emission from dust grains 
and rotational lines of carbon monoxide gas. Both methods, however, rely on unconstrained as- 
sumptions. The dust method has to assume a dust opacity per gram of dust, and grain growth can 
change this value dramatically 15 . The gas mass is then calculated by multiplying the dust mass by 
the gas-to-dust ratio, usually assumed to be ~100 from measurements of the interstellar medium^. 
The gas mass thus rests on a large and uncertain correction factor. The alternative is to use rota- 
tional CO lines as gas tracers, but these are optically thick, and therefore trace the disk surface 
temperature, as opposed to the midplane mass. The use of CO as a gas tracer then leads to large 
discrepancies between mass estimates for different models of TW Hydrae (from 5 x 1CT 4 M to 
0.06 M ), even though each matches a similar set of observations^^. 

Using the Herschel Space Observatory 17 PACS Spectrometer^ we robustly detected (9a) the 
lowest rotational transition, J = 1 — > 0, of hydrogen deuteride (HD) in the closest (D ~ 55 pc) 
and best studied circumstellar disk around TW Hydrae (Fig. 1). This star is older (3-10 Million 
yea ri 9 ! 10 !^ !) than most stars with gas-rich circumstellar diskP. The abundance of deuterium atoms 
relative to hydrogen is well characterized via atomic electronic transitions to be x-q = 1.5 ± 0.1 x 
10~ 5 in objects that reside within ~ 100 pc of the SurP^. Placing atoms into H 2 and HD, which 
is appropriate for much of the disk mass, provides an HD abundance relative to H 2 of xhd = 
3.0 x 10~ 5 . We combine the HD data with existing molecular observations to set new constraints 
on the disk mass within 100 AU - the most fundamental quantity that determines whether planets 
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can form. The disk mass also governs the primary mode of giant planet formation, either through 
core accretion or gravitational instability 21 . In this context we do not know whether our own 
Solar System formed within a typical disk, as nearly half of current disk mass estimates fall below 
the minimum mass solar nebula 8 . Our current census of extra-solar planetary systems furthermore 
suggests that even larger disk masses are necessary to form many of the planetary systems seen^ED. 

With smaller rotational energy spacings and a weak electric dipole moment, HD J = 1 — > 
is a million times more emissive than H 2 for a given gas mass at a gas temperature of T gas = 20 
K. The HD line flux (Fi) sets a lower limit to the H 2 gas mass at distance D (see Supplementary 
Information): 
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If HD is optically thick or D is hidden in other molecules such as PAHs or molecular ices, the 
conversion from deuterium mass to hydrogen mass will be higher and thus the mass will be larger, 
hence the lower limit. The strong temperature dependence arises from the fractional population 
of the J=l state which has a value of /j =1 ~ 3exp(— 128.5 K/T gas ) for T gas < 50 K in thermal 
equilibrium. Due to the low fractional population in the J=l state, HD does not emit appreciably 
from gas with T ~ 10 — 15 K, the estimated temperature in the outer disk mass reservoir (R > 
20-40 AU). The HD mass derived from Eqn. (1) provides an estimate of the mass in warm gas 
gas, and is therefore a lower limit to the total mass within 100 AU. 

The only factor in Eqn. (1) that could lower the mass estimate is a higher T gas . The upper 
limit on the J = 2 — > 1 transition of HD (Fig. 1) implies T gas < 80 K in the emitting region. This 
T gas estimate yields M gas disk > 2.2 x 10~ 4 M , but T gas is unlikely to be this high for the bulk of 
the disk. CO rotational transitions are optically thick and the level populations are in equilibrium 
with T gas , so they provide a measure of T gas - ALMA observations of CO J = 3 — > 2 emission in 
a 1.7" x 1.5" beam (radius ~ 43 AU) (see Supplementary Information and Supplementary Figure 



1) yield an average T gas = 29.7 K within 43 AU, and M gasdisk > 3.9 x 10~ 3 M . This value 
is still likely to be too low because the optically thick CO emission presumably probes material 
closer to the surface than HD, and this gas will be warmer than the HD line emitting region. Thus 
essentially all correction factors would push the mass higher than this conservative limit, which 
already rules out a portion of the low end of previous mass determinations. 

To determine the mass more accurately, we turn to detailed models that incorporate explicit 
gas thermal physics providing for substantial radial and vertical thermal structure. Both published 
models of the TW Hya disk reproduce a range of gas phase emission lines, but in one case with 
M-gas-disk = 0-06 M 14 and in the other with M gas _ disk = 0.003 M Q ^ (see Supplementary In- 
formation and Supplementary Table SI). These models were both placed into detailed radiation 
transfer simulations. The results from this calculation and the adopted physical structure are given 
in Fig. 2 for the model with M gas _ disk = 0.06 M . Fig. 2c shows the cumulative flux as a function 
of radius for the higher mass model; over 80% of the emission is predicted to arise from gas inside 
of 80 AU. Furthermore Fig. 2d provides a calculation of the HD emissive mass as a function of 
gas temperature. This calculation suggests gas with temperatures of 30-50 K is responsible for the 
majority of the HD emission. 

The model with M gas _ disk = 0.003 M predicts an HD line flux of F t = 3.8 x 1(T 19 W m~ 2 , 
more than an order of magnitude below the detected level. For this model to reach the observed flux 
the disk mass must be increased by a factor of 20, ruling out this low disk mass. The M gas - disk = 
0.06 M model predicts Fi = 3.1 x 1CT 18 W m~ 2 , still a factor of two below the observed value. 
Even the "high" mass estimate is too low. Based on this model we estimate the disk gas mass 
within 80 AU, where the majority of HD emissions arise, is 0.056 M . Both of these models 
match other observations as the low disk mass model matches CO and 13 CO J = 3 — > 2, while the 
higher mass model reproduces CO J = 2 — > 1, J = 3 — > 2, J = 6 — > 5. Both compare emission 
predictions to other species as well. However, they differ by a factor of 10 in predicting the HD 
emission. This difference shows the value of HD in constraining masses. 
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The age of TW Hya is uncertain. The canonical age of the cluster is lOl^ Million years 9 . 
However there could be an age spread in cluster members and ages estimated for TW Hya itself 
range from 3-10 Myit^E^l. Even at the low end of this range TW Hya is older than the half-life 
of gaseous disks, inferred to be about 2 Myr 8 . In the case of the TW Hya association there is 
also little evidence for an associated molecular cloud 24 , which is an additional indicator that this 
system is relatively older than most gas-rich disks. The lifetime of the gaseous disk is important 
as it sets the available timeframe for the formation of gas giants equivalent to Jupiter or Saturn. 
Based on our analysis, TW Hya contains a massive gas disk (> 0.06 M ) that is several times the 
minimum mass required to make the planets in our solar system. Thus, this "old" disk can still 
form a planetary system like our own. 

The recent detection of cold water vapor from TW Hya found indirect evidence for a large 
water ice reservoir (~ several thousand Earth oceans) assuming a disk mass of 0.02 M 25 . Our 
higher mass estimate implies a larger water ice reservoir, perhaps increased by a factor of two. The 
mass estimate in this one system lies on the upper end of previous mass measurements^, hinting 
that other disk masses are underestimated. The main uncertainty in the masses derived here is the 
gas temperature structure of the disk. Looking forward, observations of optically thick molecular 
lines, particularly CO, can be used to trace the thermal structure of gas in the disk. Observations 
of rarer CO isotopologues then provide constraints on the temperature in deeper layers 26 . With the 
Atacama Large Millimeter Array, we will readily resolve multiple gas temperature tracers inside 
80 AU where HD strongly emits. When used in tandem with HD we will be able to derive the 
gas mass with much greater accuracy (our simulations suggest to within a factor of 2-3). More- 
over, additional HD detections could be provided by Herschel and with higher spectral resolution 
by SOFIA GREAT under favorable atmospheric conditions. These data could be used alongside 
emission from species such as C 18 0, C 17 0, or the dust to calibrate these more widely available 
probes to determine the disk gas mass. Thus, with the use of HD to complement other observa- 
tions and constrain models we may finally place useful constraints on one of the most important 
quantities that governs the process of planetary birth. 
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Figure 1 : Herschel detection of Hydrogen Deuteride in the TW Hya protoplanetary disk, (a) 

The fundamental J = 1 — > line of HD lies at ~1 12/xm line. On 20 November 201 1 it was detected towards 
the TW Hya disk at the level of 9a. The total integrated flux is 6.3 ± 0.7 x 10~ 18 W m~ 2 . We also report a 
detection of the warm disk atmosphere in CO J = 23 — > 22 with a total integrated flux of 4.4 ± 0.7 x 10~ 18 
W m" 2 . The J = 1 — > line of HD was previously detected in a warm gas cloud exposed to radiation 
from nearby star^by the Infrared Space Observatory. Other transitions have also been detected in shocked 
regions associated with supernova and outflows from massive starj 28 ! 29 ! (b) Simultaneous observations of 
HD J = 2 1 are shown. For HD J = 2 ^ 1 we find a detection limit of < 8.0 x 10~ 18 W m" 2 (3a). We 



also report a detection of the OH 2 I1] 



I doublet near 55.94 /xm with an integrated flux of 4.93 ± 0.27 



x 10 W m . The spectra include the observed thermal dust continuum of ~ 3.55 Jy at both wavelengths. 
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R (AU) Gas Temperature (K) 

Figure 2: Model of the physical structure and HD emission in the TW Hya circumstellar disk. 

(a) Radial and vertical distribution of the H2 volume density, nn 2 , calculated in a model disk with mass 0.06 
MgJ^I. Contours start from the top at Logio [nH 2 /cm -3 ] = 1.0 and step downwards in units of factors of ten. 

(b) The gas temperature structure as derived by the thermochemical model 14 . Contours are 10, 25, 50, 75, 

100, 150, 200, 250, and 300 K. (c) Radial and vertical distribution of the HD J = 1 volume density, 71hd j=i> 

predicted in a model disk with the gas density and temperature structure as given in panels (a) and (b), with 

an HD abundance relative to H2 of 3.0 x 10 -5 . Contours start at the top with Logio [nnD j=i/cm~ 3 ] = —3 

and are stepped in factors of 10. The red line on the plot shows the cumulative flux contribution as a function 

of radius in terms of fractions of the overall predicted flux of 3.1 x 10~ 18 W mr 2 . To predict the HD line 

emission we calculate the solution of the equations of statistical equilibrium including the effects of line 

and dust opacity using the LIME code^. (d) Fraction of the HD emission arising from gas with different 

temperatures, computed as function of the mass of HD excited to the J=l state in gas at temperatures binned 

in units of 5 K and normalized to the total mass of HD in J=l. In particular, J tt-hd J=i2vr Rdrdz is computed 

successively in gas temperature bins of 5 K and then normalized to the total mass of HD in the J = 1 state. 
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Supplementary Information 

1. Herschel Observations 

TW Hya (a(2000) = ll h 01 m 51.91"; (5(2000) = -34 o 43'17.0") was observed as part of an open time Her- 
schel Program (OTl_ebergin_4) using the PACS 18 instrument on November 20, 2011. The PACS Range 
scan chop-nod mode was used. The background emission from the the telescope and sky was subtracted 
using two nod positions 1.5' distant from the source in each direction. The data cover a small spectral range 
(111.2-114.2 /xm and 55.6-57.1 /im), with high sampling density (for a total of 81 steps). The first range 
includes HD J = 1 — > and the second was designed to detect HD J = 2 — > 1. We used 12 repetitions to 
increase sensitivity for a deep scan, for a total integration time of 25124 seconds (36 repetitions) on TW 
Hya. The predicted line RMS was 0.62 x 10~ 18 W mT 2 and 2.65 x 10~ 18 W mT 2 at 112 and 56 ^m, 
respectively. 

The data were reduced using the Herschel Interactive Processing Environment (HIPE) 31 v8.1 pipeline (cal- 
ibration set 32). PACS is a 5x5 integral field unit spectrometer with a pixel size of 9.4" x9.4"; the source 
showed no sign of extended emission beyond the point spread function and was centered within 0.2 pixels 
(within 2") of the centerpoint in each case. We used the standard "calibration block" script to reduce the 
data for optimal signal-to-noise utilizing only the central pixel of the array, and then scaled this value to 
the spectra extracted from the central 3x3 pixels. We compared this to the PSF-corrected output from the 
pipeline, and the flux levels matched to within 10%. The PSF-corrected flux was 10% higher than the 3x3 
extraction, which indicates an overcorrection and no sign of extended emission. Thus, we increased the flux 
measured from the central-spaxel-only spectrum by 10% to match that of the 3x3 extraction, yielding line 
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fluxes, Fi = (4.4±0.7) x 1(T 18 W m" 2 for CO J = 23 -> 22 (centered at 113.509 /xm indicating a possible 
blend between CO at 113.46 /im and H 2 at 113.53 jan), and = (6.3 ± 0.7) x 10~ 18 W m^ 2 for HD 
J = 1 — > (centered at 112.086 fim). The line flux was obtained using a Gaussian fit. The continuum was 
determined by a first-order polynomial simultaneous fit in a fairly tight region around the line, avoiding the 
CO transition. The HD J = 2 — > 1 transition was not detected with a 3a upper limit of Fi < 8.0 x 10~ 18 
W m~ 2 . The lines are not resolved with A/AA ~ 1000 at 112 fim and ~ 1500 at 56 /xm. The absolute 
uncertainty on flux calibration is potentially larger, with an important factor being the overall pointing and 
presence/absence of extended flux. Because TW Hya is a point source, the observations are well centered, 
and we included a PSF correction, the uncertainty in Fi is limited to 10-20%, negligible compared to other 
uncertainties. 

From the CDMS database^, the J = 1 -> line at 2674.986 GHz has E J=1 = 128.5 K with A l0 = 
5.44 x 10~ 8 s _1 . The J = 2 -> 1 line at 5331.561 GHz has E J=2 = 384.58 K with A 2 \ = 5.16 x 10~ 7 
s . The A-coefficients have been calculated using the dipole moment of 8.56 x 10" 

2. Simple Estimate of the Gas Mass 

The mass implied by a line flux of optically thin emission from an unresolved source can be derived in two 
steps. The total number of HD molecules (A/hd) is related to the line flux by this relation, assuming that the 
beam encompasses the source: 



In this expression f u = 3.0 * exp(— 128.5 K/T)/Q(T) is the fractional of HD molecules in J = 1, 
D is the distance, and u is the frequency. Converting to mass, and assuming all is in H2, M gas disk = 
2.37 * JTIhA/hd/^HDj where xhd is the abundance of HD relative to H2, mn is the mass of a hydrogen 
atom, and 2.37 is the mean molecular weight per particle, including helium and heavy elements 34 . This 
gives 



AttD 2 



u 



(2) 




(3) 
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This estimate represents a lower limit because the HD emission may not trace all the mass in the disk. 

3. CO Emission and T gas from TW Hya 

In the main text we used the resolved ALMA Science Verification CO J = 3 — > 2 data to set a limit on 
the gas temperature. In Supplementary Figure 1 we show the integrated emission map with a beam size 
1.7" xl.5", which corresponds to a radius of 47 x 41 AU. The map demonstrates that the CO emission 
is resolved, so we assume unity filling factor. Based on the central spectrum the observed peak radiation 
temperature is Tr = 22.2 ± 0.1 K (we note that this data has a calibration uncertainty of 10%). Tr 
is linearly related to the intensity, so a correction for the difference between the Planck function and the 
Ray leigh- Jeans approximation results in an average gas kinetic temperature of T gas = 29.7 K in the beam, 
assuming optically thick and thermalized emission. Observations of CO J = 6 — > 5 confirm this value. The 
observed CO J = 6 — > 5 peak intensity is Tr = 16.9 ± 3.2 K^, which corresponds to a gas temperature of 
30.6 K. 

Assuming T gas = 30 K yields a minimum disk mass of 3.9 x 10~ 3 M . 

4. HD Emission Line Models 

For more detailed models we use two predicted physical structures (Gorti et al. and the Thi et al. TW 
Hya models) derived from thermochemical calculations that match a range of gas phase emission lines. The 
Gorti et al. TW Hya model has a disk gas mass of 0.06 M Q , and we used direct output from the model. 
For the Thi et al. TW Hya model we use a reproduction. Specifically, to reproduce this model we adopt 
the physical parameters for the dust structure as given in their 13 Table 2. The dust optical constants are 
as prescribed 3 ^ with a dust grain size distribution from 3 x 10~ 2 /im to 1 mm and a power law index of 
3.4. These inputs were placed into RADMC 37 where we verified that the model reproduced the observed 
spectral energy distribution and the dust thermal structure given in their Fig. A.2. For this calculation the 
total disk gas mass is 0.003 M and we verified that our gas density distribution matched that provided in 
their Figure A.l. We do not model the H-H2 transition; this is appropriate for HD which will emit below 
this transition region. Finally we used the gas temperature distribution given in their Fig. A. 3 as input to our 
LIME calculations for HD. As a final check we computed the predicted emission of this reproduction for 
CO and 13 CO J = 3 — > 2 to observed values 38 and found them to be in agreement to within a factor of two. 
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To predict the HD line emission from the physical models we calculate the solution of the equations of 
statistical equilibrium including the effects of line and dust opacity using the LIME code 30 . Although the 
HD line emission is fairly close to LTE, we adopt the collision rate coefficients of HD with H2 39 . The 
collision rates are insensitive to the rotational state of the H2 collision partner 40 . The disk dust optical depth 
is determined by the dust opacity coefficient k\ and the dust mass distribution within the disk. Typical 
dust mixture^K3 suggest ki%2 ^ m ~ 30 cm 2 g _1 . The excitation model explicitly explores pumping by 
continuum radiation, which is found to be negligible when compared to the effects of the dust optical depth. 

Table S 1 provides the predictions from the specific fhermochemical model including an emission calculation 
with dust opacity at 1 12 /xm (the realistic case) and without dust opacity. The Thi et al. TW Hya model with 
the lowest mass predicts an HD emission line that is too weak. To match observations within the framework 
of the Thi et al. TW Hya model we require an increase in mass by a factor of 20. We stress that this is 
approximate as scaling the mass will change the thermal solution (gas temperature, chemistry). Thus, this 
factor only illustrates the mass dependency and is not definitive. The thermochemical model that comes 
closest to the observed flux is the Gorti et al. TW Hya 0.06 M model - which is still about a factor of two 
smaller than the observed HD emission. 



Supplementary Table SI: Specific Model Predictions 



Disk Model 


Gas Mass 


HD J = 1 a 


HD J = 2 -> 1 




(M ) 


(W m" 2 ) 


(W rrr 2 ) 


Thi et al. TW Hya 


0.003 


3.8 xirr 19 


1.4 xicr 19 


Gorti et al. TW Hya 


0.06 


3.1 xicr 18 


3.3 xicr 18 


Observations 




6.3±o.7 xirr 18 


<8.o xicr 18 



Q Fluxes without dust optical depth are 7.4 x 1CT 19 W rrr 2 (Thi et al. 0.003 M ), 
4.2 x 10" 18 W rrr 2 (Gorti et al. 0.06 M ) 
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R.A. offset (arcsec) 



Supplementary Figure 3: ALMA Science Verification Observations of CO J = 3 — > 2 in TW 
Hydra. (Left) Map of CO J = 3 — > 2 integrated emission with intensity scale given on the left. 
The beam size of this observation is 1.7" x 1.5" and is shown in the figure. (Right) Blow-up 
of the central spectrum. 
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